#DESCRIPTION:
#Module computing Mohr-Coulomb shear strength parameters from a given set of Hoek-Brown strength parameters and slope height

#REFERENCE
#Hoek, E., Carranza-Torres, C., Corkum, B. (2002). Hoek-brown failure criterion - 2002 edition. NARMS-TAC (pp. 267-273).


#1. libraries and modules
import math
import numpy as np
from scipy.interpolate import UnivariateSpline
import pandas as pd
import matplotlib.pyplot as plt

#2. math/trigo functions
tan=math.tan
rad=math.radians
pow=math.pow
exp=math.exp

#3. Generating Hoek-brown failure envelope as spline defined by normal stress vs shear stress (sign vs tau)
def HoekBrownFE_spline(GSI,D,mi,sig_ci,uw,H):
    
    #1. computing Hoek-Brown derived values
    GSI=float(GSI)
    mi=float(mi)
    sig_ci=float(sig_ci)
    try:
        mb=mi*exp((GSI-100)/(28.0-(14*D)))
        s=exp((GSI-100)/(9.0-(3*D)))
        a=0.5+((pow(exp(1),(-GSI/15.0))-pow(exp(1),(-20/3.0)))/6)
        sig_c=sig_ci*pow(s,a)
        sig_t=-s*sig_ci*1.0/mb
        sig_cm=sig_ci*(mb+4*s-a*(mb-8*s))*(0.25*mb+s)**(a-1)/(2*(1+a)*(2+a))
        sig3max=0.72*sig_cm*(sig_cm/(uw*H))**-0.91
    except: 
        print "error computing derived values"

    #2. creating Hoek-Brown failure envelop as sigma3 vs sigma1
    sig3=np.concatenate((np.linspace(sig_t*.9999,-0.9999*sig_t,1000),np.linspace(-sig_t,2*sig3max,200)))
    sig1=sig3+sig_ci*(mb*sig3/sig_ci + s)**a
    ds1s3=1+a*mb*np.power(((mb*sig3/sig_ci)+s),(a-1))

    #3. deriving corresponding sign and tau values from sigma3 vs sigma1
    sign=(sig1+sig3)/2. - ((sig1-sig3)/2.) *(ds1s3-1)/(ds1s3+1)
    tau=(sig1-sig3)*np.sqrt(ds1s3)/(ds1s3+1)

    #4. Generating spline defined by sign vs tau
    sign_tau=UnivariateSpline(sign,tau,s=0)
    #plt.plot(sig3, sig1)
    #plt.plot(sign, sign_tau(sign))

    #returns spline and minimum sign value allowable
    return sign_tau, sig_t
    
#4. Computing Mohr-colomb parameters, phi and cohesion, from spline defined by sign vs tau
def phi_coh_from_HoekBrownFE_spline(sig_n_list, sign_tau,sig_t):
    phi=np.where(sig_n_list>sig_t,np.arctan(sign_tau(sig_n_list,1)),np.nan*np.ones(len(sig_n_list)))
    coh=np.where(sig_n_list>sig_t,(sign_tau(sig_n_list)-sign_tau(sig_n_list,1)*sig_n_list),np.nan*np.ones(len(sig_n_list)))
    #returns phi (in units of radians) and cohesion (in units of sig_ci in function HoekBrownFE_Spline)
    return phi,coh











    
